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Abstract 

We use first-principles techniques to re-examine the suggestion that transitions seen in high- 
P experiments on Mo are solid-solid transitions from the bcc structure to either the fee or hep 
structures. We confirm that in the harmonic approximation the free energies of fee and hep 
structures become lower than that of bcc at P > 325 GPa and T below the melting curve, as 
reported recently. However, we show that if anharmonic effects are fully included this is no longer 
true. We calculate fully anharmonic free energies of high-T crystal phases by integration of the 
thermal average stress with respect to strain as structures are deformed into each other, and also 
by thermodynamic integration from harmonic reference systems to the fully anharmonic system. 
Our finding that fee is thermodynamically less stable than bcc in the relevant high-P/high-T 
region is supported by comparing the melting curves of the two structures calculated using the 
first-principles reference-coexistence technique. We present first-principles simulations based on 
the recently proposed Z method which also support the stability of bcc over fee. 

PACS numbers: 64.10. +h,64.70.D-,64.70.K-,71.15.Pd 
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I. INTRODUCTION 



The past 10 years have seen a lively controversy over the phase diagrams of transition 
metals at megabar pressures, with the pressure dependence of the melting temperature 
dT m /dP from diamond-anvil-cell (DAC) measurements differing greatly from that deduced 
from shock data and from first-principles calculations (see Fig. [T]). lr — One suggested reso- 
lution of the controversy is that the transition interpreted as melting in some of the DAC 
experiments may in fact be a solid-solid crystallographic transformation, and in the case 
of molybdenum this is consistent with the observation of two transitions in shock experi- 
ments.— 1 ^ This suggestion appeared to be confirmed by recent first-principles work on Mo, 
which indicated a transition from the low-temperature body-centred-cubic (bcc) structure 
to a close-packed structure in the appropriate temperature region.- 1 ^ 1 ^ However, that work 
relied on two important assumptions, which we examine in detail in this paper. The results 
we shall present imply that those assumptions and the conclusions drawn from them may 
be incorrect, so that further work is still needed to resolve the controversy. 

DAC measurements have been reported on the melting curves of several transition metals, 
including Ti, V, Cr, Fe, Co, Ni, Mo, Ta and W.—&^£- The measurements extend up to nearly 
100 GPa (1 Mbar), and in most cases the increase of melting temperature T m between 
ambient pressure and 100 GPa is surprisingly small; in the case of Mo, the increase is 
only ~ 200 K.-£ These findings are in stark contrast to the melting curves deduced from 
shock measurements, which are available for Fe, Mo, Ta and W.— ^ For Mo, the increase 
of T m between ambient and 100 GPa estimated from shock data is ~ 2000 K. Recently, 
density functional theory (DFT) calculations of the melting curves of Mo and Ta have been 
reported.-^ 1 ^^ The DFT predictions are expected to be reliable, because it is well known 
that DFT, without any adjustable parameters, gives excellent results for a wide range of 
properties of transition metals, including cold compression curves up to ~ 300 GPa,->2*22>2i 
Hugoniot P(V) curves,- 1 ^ 1 ^ phonon dispersion relations (and their pressure dependence, in 
the case of Fe),- 1 ^^— and low-temperature phase boundaries.— Furthermore, techniques for 
calculating melting curves using DFT have become firmly established over the past 10 years 
and more, and are known to give accurate results.— ~— The DFT results for T m (P) of Mo 
and Ta lend support to the correctness of the melting curves deduced from shock data.-^i^ 

It has become clear very recently that experimental difficulties may have led to a sub- 
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stantial underestimate of high-P melting temperatures in earlier DAC measurements. The 
work of Dewaele et al.— indicates that formation of metal carbide by chemical reaction be- 
tween the diamond and the metal sample can be a major problem. Their work also shows 
that difficulties in the pyrometric measurement of temperature can also lead to substantial 
underestimates of T m . In the case of Ta, their measurements give a melting curve that is well 
above those given by earlier DAC experiments and is fairly close to (though systematically 
lower than) the predictions from DFT. Nevertheless, in the case of Mo, the occurrence of 
two breaks in the shock data seems to leave little doubt that there is a transition from the 
low-T bcc structure to an unidentified high-T crystal structure, followed by the melting of 
the latter. The T and P of the lower transition (3500 K and 200 GPa) lie close to the natu- 
ral extrapolation of the T(P) boundary identified as melting in the older DAC experiments. 
This suggests that this boundary is associated with a bcc-solid transition. Even in the case 
of Ta, recent evidence based on DFT calculations 9 indicates that a bcc-solid transition is 
implicated in earlier DAC attempts to detect high-P melting. 

To substantiate the picture of a bcc-solid transition followed by a melting transition, it 
is necessary to show that another crystal structure becomes thermodynamically more stable 
than bcc at high temperatures, and to identify this structure. This was the aim of the 
recent DFT work on Mo by Belonoshko et al 7 They showed that, in the quasiharmonic 
approximation, the Gibbs free energy of the fee structure is lower than that of bcc over a 
substantial high-P/high-T region of the phase diagram below the bcc melting curve. The 
fee structure becomes harmonically unstable (there are imaginary phonon frequencies) for 
P < 350 GPa, but extrapolation of the predicted bcc-fcc phase boundary passes quite 
close to the (P, T) of the lower shock transition. The quasiharmonic calculations on the 
bcc and fee free energies were independently confirmed by two subsequent papers.—^ As 
independent evidence that fee is more stable than bcc at high T, Belonoshko et air- used the Z 
method^ 1 ^ to calculate the melting curve of fee Mo. (The Z method employs observations 
of spontaneous melting of the superheated solid in constant-energy molecular dynamics 
simulations to determine points on the melting curve.) They found that the fee melting 
curve lies above the bcc melting curve, thus appearing to confirm that the free energy of fee 
is lower than bcc. However, we note the two important assumptions made here: first, that 
anharmonic contributions to the free energies of high-T bcc and fee Mo can be neglected; 
and second, that the first-principles statistical-mechanical techniques that were employed 
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have the precision needed to distinguish between the possibly rather similar melting curves 
of bcc and other stuctures. There is also the question of whether fee can remain vibrationally 
stable at high T in the region P < 350 GPa, where the harmonic phonons are unstable. 
These are the issues addressed in the present paper. 

We use two methods here for comparing the free energies of the bcc and other crystal 
structures, and both methods fully include anharmonicity. The other structures we examine 
are fee and hexagonal-close-packed (hep). The first method uses the fact that the free 
energy difference between two systems that differ only by a finite strain can be obtained by 
integrating the thermal average stress with respect to strain.—^ We use this idea to obtain 
the free energy difference between fee and bcc, which can be transformed into each other by 
a continuous strain along the Bain path (BP). The second method employs thermodynamic 
integration (TI) from a harmonic reference system to the fully anharmonic system described 
by DFT. This is essentially the same method as we employed in earlier work on Fe.-*22^ 
As further ways of probing the possible thermodynamic stability of the fee structure, we 
have re-examined the melting curve of fee Mo, using the reference coexistence technique 
employed in some of our earlier work,— MdMM. an d we have also performed our own first- 
principles Z-method calculations on the melting of bcc and fee Mo. All the results point to 
the conclusion that none of the other structures is thermodynamically more stable than bcc 
at high T. 

The rest of the paper is organised as follows. In Sec. |TTJ we summarise briefly the technical 
details of the DFT methods employed in all the calculations, and we then present our results 
for the harmonic dispersion relations of the bcc, fee and hep structures over a wide range of P; 
we note the pressure thresholds below which the fee and hep structures become harmonically 
unstable. In the same Section, we report our results for the harmonic free energies, and hence 
the predicted phase boundaries separating the different structures. Sec. [Ill] presents our 
results on the vibrational, elastic and thermodynamic stability of the different structures, 
including our Bain-path calculations of the free energy differences between fee and bcc. 
Our calculations of the anharmonic contributions to the free energy, and the effect of these 
contributions on the phase boundaries are presented in Sec. IIV1 Our reference coexistence 
calculations of the fee melting curve and the comparison with the bcc melting curve obtained 
already by the same technique are outlined in Sec. |V] where we also present our new Z- 
method calculations. Finally, we draw all the results together in Sec. ED and suggest what 
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FIG. 1: Points on the solid-liquid boundary of Mo as observed in DAC (A |3l) and shock wave 
(□ [13|]) experiments. The shock wave datum (• [13]) obtained at low-P has been interpreted as a 



solid-solid phase transition, 
by Cazorla et al. (solid line 
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l[) and Belonoshko et al. (dashed line [7|). 



future investigations might help to resolve the controversies over the phase diagram of Mo 
and other transition metals. 



II. HARMONIC CALCULATIONS 



A. DFT techniques 

All calculations were done using the projector augmented wave version of DFT as imple- 
mented in the VASP package.— 1 ^ All atomic states up to and including 4s were treated as 
core states, with 4p and all higher states being valence states. We used the PBE form of gen- 
eralised gradient approximation to the exchange-correlation functional.— An energy cut-off 
of 224.6 eV was used throughout; the adequacy of this value was shown in a previous work 
where we performed extensive numerical convergence tests.- Dense Monkhorst-Pack grids^ 
were used for electronic Appoint sampling in static perfect lattice calculations, to guarantee 
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FIG. 2: Ab initio vibrational phonon frecuencies of Mo in the fee structure calculated at volumes 
V = 9.64 A 3 /atom (P = 328 GPa, solid line) and V = 10.19 A 3 /atom (P = 265 GPa, dashed line). 

convergence of the total energy to better than 1 meV/atom. Thermal excitation of electrons 
was included via the finite-T version of DFT originally developed by Mermin.— Phonon 
frequencies in our calculations were obtained by the small-displacement method^ 1 ^ using 
large supercells. For molecular dynamics (MD) simulations, we used the Born-Oppenheimer 
scheme where the self-consistent ground state is recalculated at each MD time-step. These 
simulations were performed in the microcanonical (N, V, E) and canonical (N, V, T) ensem- 
bles; temperatures in (N, V, T) simulations were maintained using Nose thermostats. Values 
of the technical parameters (duration of the MD runs, k-point grids, number of atoms, etc.) 
will be presented together with the results. 

B. Harmonic free energies and phase boundaries 

It is convenient to represent the total Helmholtz free energy F(V,T) of the system at 
volume V and temperature T as the sum of three parts:— 



F(V, T) = F P (V, T) + F h (V, T) + F^V, T) . 



(1) 



6 



N 

X 



o 

CD 

cr 
o 



CO 
CD 



i 

CD 

E 
o 

CD 

o 



10 - 



9 - 



8 - 



7 - 



BCC 


1 1 

• 




"FCC 


A 


— /m 


HCP 


▲ 






























1 1 1 



14 



13 



12 



11 

V (A 3 /atom) 



10 



FIG. 3: Geometric mean frequency ui of Mo in the bcc, fee and hep structures as a function of 
volume. Symbols represent states at which the calculations have been carried out and the lines are 
guides to the eye. 

Here, F p is the Helmholtz free energy of the static perfect lattice: it is a free energy, because 
we include thermal electronic excitations.— 1 ^ The second term is the free energy due to 
lattice vibrations, calculated in the harmonic approximation. The remainder F a accounts 
for anharmonicity; we ignore F a here, but show how to compute it in Sec. IIVI 

The calculation of F p (V,T) is completely standard. It is known from previous work 
that at T = K and pressures P < 660 GPa the most stable phase of Mo is bcc.- 1 ^ At 
higher compressions, Mo stabilizes in the double hexagonal closed-packed (dhep) structure 
as recently shown by Belonoshko et al. (see Fig. 1 of Ref. jzj]) and confirmed in our own 
calculations. Since we focus here on pressures P < 600 GPa, we have computed F p (V,T) 
on a grid of (V, T) points for the bcc, fee and hep structures. This grid spans the ranges 
8.25 < V < 15.55 A 3 /atom and < T < 10000 K, with state points taken at intervals of 
0.5 A 3 /atom and 500 K, respectively. We then fit the F P (V,T) results obtained at fixed T 
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FIG. 4: Solid-solid phase boundaries in Mo at high P and high T as obtained with first-principles 
harmonic free-energy calculations. Results obtained by Belonoshko et al. |7| are shown for com- 
parison. 



to a third-order Birch-Murnaghan equation 50 of the form 
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F P (V,T) = E + -V K 
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where E Q and K Q = —Vod 2 E/dV 2 are the values of the energy and the bulk modulus at 
equilibrium volume V , respectively, x — f (4 — -^o) an< ^ -^o = [dK/dP], with derivatives 
evaluated at zero pressure. Finally, the dependence of parameters E , K , V$ and K' on T 
is fitted to 4-th order polynomial expressions. 

We obtain the harmonic phonon frequencies u^ s by diagonalising the dynamical matrix, 
which is the spatial Fourier transform of the force-constant matrix. Our calculations of the 
latter by the small-displacement method,— 1 ^ used large supercells of 216 atoms (4x4x4 
/c-point grid) for the bcc and fee structures and 200 atoms for hep (4x4x3 fc-point grid). We 
performed extensive tests for Mo in the bcc structure which showed that these parameters^ 
guarantee F^ values converged to less than 1 meV/atom; these parameters are assumed to be 
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equally adequate for the fee and hep structures. In principle, the force-constant matrix and 
the frequencies u^ s depend on the electronic temperature, but we ignore this dependence 
here. Phonon calculations performed with an electronic T equal to 2000 and 5000 K provide 
Fh results that agree within 1 meV/atom, so we used T = 2000 K in all the cu q)S calculations. 

The phonons are stable for bee over the entire range < P < 600 GPa, as is known 
from previous work 7,16 . However, the phonons for fee and hep are stable only above a 
threshold pressure P t h- To illustrate this, we show in Fig. [2] the fee phonon frequencies 
at V = 9.64 A 3 /atom (P = 328 GPa) and V = 10.19 A 3 /atom (P = 265 GPa). We 
see that the phonon instability first occurs at a finite wavevector, which we estimate as 
q inst = (27r/a ) (1/4, 1/4,0). We find threshold pressures of P th = 310 and 325 GPa for the 
fee (V = 9.78 A 3 /atom) and hep (V = 9.69 A 3 /atom) structures. 

When calculating the harmonic vibrational free energy F^(V,T), we use the classical 
expression: 

F h (V, T) = 3k B T ln(Tiw/k B T) , (3) 
where u is the geometric mean frequency, defined by: 

iV q -^lnK,/a;) = l, (4) 

with the sum going over all iV qjS phonon modes (wavevector q, branch s) in the first Bril- 
louin zone. This classical formula for is valid at temperatures well above the Debye 
temperature, which for Mo is around 400 K at equilibrium. We have checked that even at 
T = 1000 K the difference between the Fh values obtained with the classical and quantum 
formulas is less than 1 meV/atom, so our choice does not affect the accuracy of the results. 

We show the calculated mean frequencies Co for all the structures in Fig. [3j Comparison 
of the Cj values indicates that harmonic contributions to the free-energy tend to stabilize the 
fee and hep structures over bec, since u hcp < u hcc and u fcc < u hcc in the V range studied. 
The stabilisation is slightly greater for hep than for fee. 

We fit the dependence of the quantity ln(ftw) on V to a 3rd-order polynomial for all 
the structures in order to know the value of Fh at any (V, T) thermodynamic state using 
formula ((3]). From the harmonic free energies F'(V,T) = F P (V,T) + Fh(V, T), we have de- 
termined the transition bcc-fcc and bec-hep pressures at each temperature using the double- 
tangent construction. The phase boundaries given by these calculations are shown in Fig. HJ 
where we also indicate the harmonic phase boundaries from the calculations by Belonoshko 
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et al.— (The boundaries given by the very recent quasiharmonic calculations of Zeng et a/. 16 
are similar.) In fact, there are some discrepancies. For instance, those calculated by Be- 
lonoshko et al. lie at somewhat lower T than ours, and the slopes of the two harmonic 
bcc-hcp boundaries differ in sign. Despite these differences, we agree with Belonoshko et al. 
that in the harmonic approximation the stable high-P/high-T structure of Mo is hep. 

III. VIBRATIONAL, ELASTIC AND THERMODYNAMIC STABILITY 

Do the phase boundaries predicted by harmonic theory have anything to do with the 
transitions seen in DAC and shock experiments? If they do, then the simplest hypothesis 
is that these transitions lie on a continuation of the predicted boundaries. But in order for 
this to be true, several conditions must be satisfied. First, since the lower shock transition 
occurs at P = 220 GPa, which is far below the harmonic stability limit for both fee and hep, 
the system must somehow be vibrationally stabilised, presumably by anharmonic effects. 
Second, the system must remain elastically stable at pressures P < 220 GPa, i.e. small, 
arbitrary volume-conserving strains must not cause the free energy to decrease. Third, 
the crystal structures must have lower free energies than bcc. Vibrational stability can be 
tested by straightforward first-principles MD simulations, as has been shown in our earlier 
work on high-P/high-T bcc Fe,— and in recent work by Asker et al. on low-P fee Mo;— 
we report tests here for fee Mo. The strain dependence of free energy can be also probed 
by MD calculations in which the thermal average stress is monitored. For the case of fee, 
calculation of stress as a function of strain along the Bain path also allows us to test its 
thermodynamic stability. 

A. Vibrational stability 

When we say that a crystal in thermal equilibrium is vibrationally stable, we mean that 
the thermal average position of each atom remains centred on its pefect-lattice site, and 
does not acquire a permanent deviation away from that site. To test this, it is convenient 
to use the so-called position correlation function p(t), defined by:— 

p(t) = <(r,(t + to) -Hft-Cr^-Hj)), (5) 
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FIG. 5: Calculated mean squared displacement Ar(t) 2 (dashed line) and position correlation func- 
tion p{t) (solid line) of Mo in the fee structure as a function of time. The simulations were performed 
at two different temperatures and fixed volume V = 10.50 A 3 /atom. The value of functions Ar(t) 2 
and p(t) is in units of A 2 . 
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where rj(t) is the position of atom i at time t, R° is the perfect-lattice position of the atom, 
to is an arbitrary time origin, and ( ■ ) denotes the thermal average. In practice, the thermal 
average is performed by averaging over i an d over atoms. At t = 0, p(t) is simply the 
vibrational mean square displacement. The crystal is vibrationally stable if p(t — > oo) = 0, 
because the vibrational displacements at widely separated times become uncorrelated. But 
if atoms acquire a permanent vibrational displacement, then p(t — > oo) becomes non-zero. 
The characteristic behaviour of p(t) in a vibrationally unstable crystal can be seen in our 
earlier work on high-P/high-T Fe in the bcc structure.— For this test to work, the atoms 
must not diffuse from one site to another, and we routinely test for lack of diffusion by 
monitoring the time-dependent mean square displacement Ar(t) 2 = ( | (t + t ) — rj(t )| 2 ), 
which, in the absence of diffusion, goes to a constant equal to twice the vibrational mean 
square displacement in the limit t — > oo. 

We have performed a set of first-principles MD simulations on Mo in the bcc and fee 
structures at a wide range of thermodynamic states. A typical MD run consisted of 10 4 
steps performed with a time-step of 1 fs, with the first 5 ps allowed for equilibration, and 
only the last 5 ps used to accumulate statistical averages. The simulation cell contained 125 
particles for both fee and bcc structures and T-point electronic fc-point sampling was used. 

The MD runs were carried out for a total of 20 state points, spanning the ranges 200 < 
P < 600 GPa and 1500 < T < 10500 K. In Fig. [5], we show the mean squared displacement 
Ar(t) 2 and the position correlation function p(t) calculated for fee Mo at volume V = 
10.50 A 3 /atom and T = 1500 and 4000 K. At the lower T, the system is vibrationally 
unstable, as shown by the long-t behaviour of p{t). Clearly, atomic liquid-like diffusion does 
not occur as shown by the fluctuation of Ar(t) 2 about a constant value at long t. Since these 
MD simulations are very demanding, we did not attempt to estimate an accurate boundary 
in the P — T plane separating stable and unstable states of the fee structure. Nevertheless, 
we can say that vibrational instability was not observed in our simulations at temperatures 
T > 3000 K and pressures below P t h = 310 GPa. The recent work of Asker et al— using 
techniques similar to those used here, showed that even at P ~ GPa fee Mo is vibrationally 
stable for T > 3000 K. 
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FIG. 6: Ab initio free energy calculations performed along the Bain path at zero temperature. 
Energy differences are represented with symbols, and lines are guides to the eye. 

B. Elastic and thermodynamic stability 

The elements of the stress tensor a a p in a thermal-equilibrium system can be defined 

aS (7 (yR = V- 1 (dF /d 

£<xp)ti WIiere F is the Helmholtz free energy, e Q( g is the strain tensor, 
and V is the volume. This relation can be integrated to obtain the difference of free energy 
between two states that differ by a finite homogeneous strain.—^ The bcc and fee structures 
can be continuously deformed into one another by such a strain, following the Bain path. 
This means that the free energy difference between the two structures can be obtained by 
performing a series of MD simulations along the Bain path, calculating the thermal average 
a a p in each simulation, and then integrating numerically with respect to e a p. A necessary 
condition for elastic stability of the fee phase is that F must be a local minimum along the 
Bain path.— 

The Bain path is based on the idea that the bcc and fee structures can be regarded as 
special cases of the body-centered tetragonal lattice (bet, IA/mmm space group). Taking 
primitive vectors &i = (1, 0, 0)a, a 2 = (0, 1, 0)a, a 3 = (1/2, 1/2, ()a, the values ( = 1/2 and 
and C = 1/a/2 correspond to bcc and fee, respectively. By varying Cfrom 1/2 to l/y/2, while 
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FIG. 7: Top: Stress tensor components calculated at different C-points along the Bain path at V = 
8.26 A 3 /atom and T = 9000 K. Statistical errors are represented with bars equivalent to 0.4 GPa. 
bottom: Free energy difference AF(C) obtained at (V, T) states (8.26, 9000) = •, (10.08, 6000) = A, 
(10.54,4000) = o and (11.00,4000) = ▲ in units of A 3 /atom and K, respectively. Numerical 
uncertainties are represented with bars of 5 meV/atom and the lines are guides to the eye. 
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varying a so as to keep the volume ct 3 £ of the unit cell constant, one structure is transformed 
continuously into the other. If we denote by i*bct(C) the free energy for a given £ value, then 
the work done on going from the bcc value £ = 1/2 to the another value at constant volume 
is readily shown to be: 

F bct {() - F hcc = \v [ C [a xx + a yy - 2a zz )^d(' . (6) 

6 Jl/2 (, 

For ( = 1/V2, we obtain the free energy difference of interest F fcc — F hcc . 

As a preliminary test of the correctness of our procedures, we have performed calculations 
at T = K, in which case the free energy difference at constant volume is simply the energy 
difference. We show in Fig. [6] the results of integrating the stress for a range of ( values, 
starting from bcc. As expected, at P = 550 GPa the difference AE = Ef cc — E\y CC has the 
small value 0.068 meV/atom; at P = 350 GPa, which is close to the pressure at which the 
fee structure becomes elastically stable, AE has the much larger value 0.354 meV/atom, 
and the slope of AE is close to zero at the fee structure; at P = GPa, the curvature of 
AE is downwards, so that the fee structure is elastically unstable. 

Before starting full DFT Bain-path calculations, we have made preparatory tests to find 
out how to design the simulations so as to obtain useful accuracy. These tests were done with 
an embedded-atom empirical potential (EAM),— ^ which was tuned to reproduce the ener- 
getics of Mo in the bcc and fee structures as described by DFT MD simulations performed at 
V = 9.64 A 3 /atom and T = 7500 K. The values of the corresponding EAM parameters are, 
with the same notation as in Ref. [l| (Eq. (1)), e = 0.2218 eV, a = 5.5525 A, C = 4.3164, 
n = 3.33 and m = 4.68. We set the requirement that integration along the Bain path should 
give the free energy difference AF = Ff cc — Fb cc with errors of no more than ~ 10 meV 
due to statistical uncertainty, number of ^-points for numerical integration, and system size. 
Our tests indicated that the statistical uncertainty in o a p should be less than ~ 0.5 GPa, 
and this is achieved with runs of 3 — 4 ps after equilibration. Using the trapezoidal rule 
for numerical integration, we find that nine ( points (including the end-points ( = 1/2 and 
1/V2) suffice. 

Guided by the results of these tests, we performed the DFT MD calculations on systems 
of 125 atoms (r-point sampling), at nine ( values, with an equilibration time of 2 ps and a 
statistical sampling time of 3 — 4 ps; we use our standard plane-wave cut-off of 224.6 eV. 
(Checks on the adequacy of r-point sampling and our standard cut-off are noted below.) The 
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Bain-path calculations were done at four different (V, T) states (units of A 3 /atom and K): 
(8.26, 9000), (10.08, 6000), (10.54, 4000) and (11.00, 4000), where the pressures corresponding 
to the bcc structure are 600, 283, 226 and 187 GPa, respectively. As an illustration, we show 
in Fig. 0the computed values of a xx , a yy and a zz as a function of ( for the state (8.26, 9000). 
The Figure also reports the free energy difference Fb c t (C) — -^bcc obtained by integration at 
these four states. 

Two important conclusions are clear from these results. First, the fee structure is ther- 
modynamically unstable with respect to bcc at all the high-T states we have examined. This 
is true even at the state (600 GPa, 9000 K), which lies on the bcc- fee boundary predicted by 
harmonic theory, and at the state (283 GPa, 6000 K), which is somewhat above the exten- 
sion of that boundary. The second conclusion is that fee appears to be elastically unstable 
for the three states having P < 300 GPa, though it is weakly stable at (600 GPa, 9000 K). 
The conclusions appear to be robust, since they would be unchanged even if the statistical 
errors were considerably greater than those we have achieved. We have tested for the pos- 
sible effect of systematic errors coming from our use of T-point sampling and our standard 
plane- wave cut-off. To test the effect of cut-off, we have repeated the calculation of the ther- 
mally averaged stress components a xx , a yy and a zz at £ = 0.60 for the thermodynamic state 
V = 8.26 A 3 /atom and T = 9000 K, using the increased plane-wave cut-off of 280.7 eV. We 
find that this increased cut-off leaves all the stress components unchanged within ~ 1 GPa. 
We have done the same thing with the standard cut-off but now using Monkhorst-Pack 
(2x2x2) sampling of 8 fc-points. This has the effect of shifting all three stress components 
down by ~ 4 GPa. Since they are all shifted by essentially the same amount, this does not 
affect the integral of eqn (jBJ), so that the free-energy difference between fee and bcc remains 
unchanged.— 

Since the Bain-path calculations fully include anharmonicity, and since they are com- 
pletely at odds with the harmonic predictions, it appears that anharmonic contributions to 
the free energy must be very substantial at high temperatures. We examine these contribu- 
tions directly in the next Section. 
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FIG. 8: Averaged (Ui — Uq)\ values obtained at different A-points in anharmonic free energy 
calculations perfomed for bcc Mo at V = 10.08 A 3 /atom and T = 6000 K. The dashed line 
corresponds to a 4-th order polynomial curve used to reproduce the variation of these values on 
parameter A. 

IV. ANHARMONIC FREE ENERGY 

We calculate the anharmonic contribution to the free energy using thermodynamic inte- 
gration, which we have used extensively in previous work on the free energy of transition 
metals.- 1 ^ 1 ^ The general principle is that we compute the change of Helmholtz free energy 
as the total energy function U\(r 1 , . . .r N ) is changed continuously from U to U\, the free 
energies associated with these energy functions being F and F\. Then the thermodynamic 
integration formula is: 

F X -F Q = [ (U 1 -U )xd\, (7) 
Jo 

where (-)x is the thermal average evaluated for the system governed by the energy function 
U\ = (1 — \)Uq + XUi. In practice, we take Ui to be the DFT total energy function U, whose 
free energy we wish to calculate, and Uq to be the total energy function U re f of a "reference" 
system, chosen so that its free energy F re f can be evaluated exactly. Here, we choose the 
reference system to be a perfectly harmonic system^. For a volume where the harmonic 
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FIG. 9: Dependence on temperature of the free energy difference between bcc and fee Mo as 
given by anharmonic and quasi-harmonic DFT free energy calculations performed at two different 
volumes. 

phonons are all stable, we can choose C/ re f to be the total energy of the DFT system calculated 
in the harmonic approximation. For V where DFT gives imaginary phonon frequencies, 
the total free energy cannot be separated into perfect-lattice, harmonic and anharmonic 
components (Eq. (JTJ). However, we know from Sec. IIII Al that the fee (and hep) system can 
still be vibrationally stable at such volumes, so that the free energy should still be calculable. 
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In these cases, we create a harmonic reference system by adding artificial on-site harmonic 
springs to remove the harmonic instability. 

In applying this scheme in practical DFT calculations, there is a subtle point connected 
with electronic fc-point sampling, which we note here. Ideally, we should use infinitely fine 
/c-point sampling; we denote the DFT total energy calculated in this way by U°°(ti, . . . rjy). 
(As usual, U°° is a free energy, because it includes thermal electronic excitations.) However 
practical DFT simulations have to be performed with limited /c-point sampling, and we 
denote the total energy in this case by U k {ji ) . . . tn). (In fact, most of our simulations are 
performed with T-point sampling.) Now with both perfect and imperfect fc-point sampling, 
we can separate the total energy into the total (free) energy of the perfect lattice U v and the 
vibrational energy f/ vib . We write U°° = + U™ h and U k = U k + U k ih . Now the energies 
U£° and are very large, and we do not wish to incur fc-point errors in these; we do not 
need to do so, since £/£° and U k can be calculated explicitly in advance. In a practical DFT 
simulation with limited Appoint sampling, we therefore make smaller errors if we take the 
total energy to be U£° + U k ih = U£° + {U k — U k ). The reference system should be taken to 
have the total energy U ref = U£° + U^ , where t/£ is a bilinear function of the displacements 
of atoms from their regular lattice sites. 

With these points in mind, the A-dependent total energy function used in thermodynamic 
integration is: 

U x = C/ p °° + (1 - W f + \(U k - U k ) . (8) 
The total free energy of the system is then: 

F = U™ + + f dX (U k - U k - Uf) x . (9) 
Jo 

We evaluate the integral in Eq. (J2J) numerically. For this, first we perform a series of 
ab initio molecular dynamics simulations in the (N, V, T) ensemble governed by the energy 
function U\ at different A-values. We then fit a fourth-order polynomial to the (U k — U k — 
U^)\ values obtained from these simulations, in order to perform the A-integration. Our 
tests show that DFT MD simulations performed at five equally spaced A-points are enough to 
ensure convergence of the free energy to better than 10 meV/atom (see Figure E}. A typical 
run consisted of 3 x 10 3 MD steps performed with an time-step of 1 fs with statistical averages 
taken over the last 2 ps. This procedure gives values for (U k — U k — U[ ci )\ converged to 
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V T 


pbcc 


pbcc 


pfcc 


F fcc 


phcp phcp 


9.20 1000 


-5.694 


-0.002 


-5.447 


-0.005 




(390 < P < 415) 3000 


-7.014 


-0.005 


-6.825 


-0.004 


-6.619 -0.015 


4500 


-8.295 


0.006 


-8.156 


0.020 


-7.964 0.003 


6000 


-9.747 


0.001 


-9.641 


0.062 


-9.472 0.022 


7500 


-11.336 


-0.017 


-11.258 


0.104 


unstable 


9000 


-13.047 


-0.051 


-12.979 


0.149 




9.64 3000 


-8.117 


0.000 


-7.885 


-0.008 




(325 < P < 345) 4500 


-9.409 


-0.005 


-9.238 


0.035 




6000 


-10.872 


-0.005 


-10.742 


0.098 




7500 


-12.484 


-0.038 


-12.367 


0.170 




10.08 3000 


-9.031 


0.005 


-8.783 






(270 < P < 285) 4500 


-10.338 


0.001 


-10.160 






6000 


-11.819 


-0.019 


-11.691 







TABLE I: Total and anharmonic free energy values obtained for Mo in different crystal structures 
within the thermodynamic range 270 < P < 415 GPa and 3000 < T < 9000 K. The cases in which 
the value of the F a = F — F p — Fh term is not shown correspond to thermodynamic states at which 
the corresponding crystal structure is harmonically unstable. Volumes are in units of A 3 /atom, 
pressures of GPa and free energies of eV/atom. 

better than 6 meV/atom. The simulation box employed contains 125 particles (128 in the 
hep case) and T-point electronic sampling was used. 

Anharmonic free energy results are shown in Table I. We see that solid Mo is always 
thermodynamically more stable in the bec structure than in the other structures examined. 
This conclusion disagrees with DFT calculations performed in the quasi-harmonic approx- 
imation (Sec. [II]) which predict fee and hep Mo as more stable than bec Mo at pressures 
and temperatures above ~ 350 GPa and ~ 5000 K. Moreover, the total free energy of Mo 
in the hep phase is around ~ 0.1 eV/atom larger than in the bec or fee structures. The 
reason behind these results is that the anharmonic energy term F a in general is negative 
for the bec structure while positive for the rest of structures, particularly at low pressures 
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FIG. 10: Ab initio (Al) high-P and high-T melting curve of Mo calculated for the bcc [1[ (Cazorla 
et al.) and fee (present work) crystal structures. Melting (P™ , ) states obtained in the two- 
phase coexistence simulations performed with EAM potentials are also displayed. A [J and □ Q 
represent DAC and shock-wave data, respectively. The melting line of fee Mo as calculated by 
Belonoshko et al. [7( is shown for comparison. 

and high temperatures (see Table I and Fig. |H]). We find very good agreement between the 
results obtained using thermodynamic integration (TI) and integration of the stress tensor 
with respect to strain along the Bain path (BP); for instance, at state (10.08, 6000) the free 
energy difference Fb cc — Ff cc obtained with TI is —0.128 eV/atom while the corresponding 
BP value is —0.117 eV/atom. 

The results of the present Section and the previous one all indicate that neither fee nor 
hep becomes thermodynamically more stable than bcc at high temperature. If this is true, 
then the melting curves of those two crystal structures should lie lower in temperature than 
the bcc curve. We turn to this question for fee Mo in the next Section. 
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V. MELTING CURVE OF FCC MO 



There are several well established techniques for calculating first-principles melting 
curves,— including the calculation of the free energies of solid and liquid using thermo- 
dynamic integration from reference systems, and the direct first-principles simulation of 
coexisting solid and liquid in large systems. Here, we begin (Sec. IV Aj) by using the "refer- 
ence coexistence" method,-^ 1 ^^ 1 ^, because it is fairly easy to apply and because we have 
used it recently to determine the DFT melting curve of bcc Mo. We shall see that the results 
given by this method are inconsistent with earlier results for the relation between the melt- 
ing curves of bcc and fee Mo obtained by the Z method. 7 In order to investigate the reasons 
for this discrepancy, we shall present (Sec. IV Bp our own DFT Z-method calculations, using 
larger simulated systems and longer simulation times than were used in the earlier work. 

A. Reference coexistence 

The reference coexistence technique consists of three steps. First, an empirical reference 
model is fitted to first-principles simulations of solid and liquid at thermodynamic condi- 
tions close to the expected melting curve. Next, the reference model is used to perform 
simulations of coexisting solid and liquid in large systems consisting of many thousands of 
atoms, so as to find points (P™ f ,T™ f ) on the reference melting curve. Finally, differences 
between first-principles and reference free energies of the solid and liquid are used to esti- 
mate the differences between reference and first-principles melting curves. In the case of 
Fe, reference coexistence results have been compared with melting curves obtained both by 
first-principles free energy calculations and by direct first-principles simulation of coexisting 
solid and liquid, and the agreement was excellent.- 1 ^ Moreover, notable agreement between 
reference coexistence results and diffusion Monte Carlo free energy calculations has been 
also proved recently.— 

The reference model used in our reference coexistence calculations on the melting of bcc 
Mo was an embedded atom model (EAM), details of which are given in Ref. jy. We use 
exactly the same model with the same parameters here. In our work on bcc Mo, we showed 
that EAM coexistence simulations on cells containing 6750 atoms give accurate reference 
melting curves, and we use the same size of system here. The protocols used to prepare 
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the two-phase system are the same as those used before, and we accept a thermodynamic 
state {P^, T™ f ) as lying on the reference melting curve if the two phases remain in stable 
coexistence for 50 ps or more. The reference melting curve obtained for fee Mo in the present 
work is compared with our published reference curve for bec Mo in Fig. [TU1 The two curves 
are essentially identical. 

The leading-order shift AT m in melting temperature caused by going from the reference 
to the first-principles total-energy function is: 

AT m = AG ls (T£ { )/Sl s e{ . (10) 

Here, AG ls = AG 1 — AG S , where AG' and AG S are the isobaric-isothermal changes of 
Gibbs free energy of liquid and solid due to the change AU of total-energy function; the 
denominator S^ f is the reference entropy of fusion, i.e. the difference between the entropies 
of liquid and solid in the reference model. The free energy shifts AG 1 and AG S are calculated 
using the formula: 

AG = (AU) Iei - ^/3(6AU 2 ) rc{ - 1 -Vk t AP 2 , (11) 

with (3 = 1/ksT, 5AU = AU — (AU) Te { (averages taken in the reference system), kt is the 
isothermal compressibility and AP is the isochoric-isothermal difference of pressure between 
first-principles and reference systems. 

Following the procedures used in our work on bec Mo, we evaluated S^ f and the reference 
Kt values for solid and liquid using separate solid- and liquid-state simulations on cells of 
3375 atoms at (P, T) points on the reference melting curve. The values of (AZ7) re f , (5AU 2 ) TC { 
and AP were obtained from solid- and liquid-state simulations on systems of 125 atoms, 
using a 2 x 2 x 2 Monkhorst-Pack grid for electronic /c-point sampling. 

In Fig. dni we compare the resulting DFT melting curve for fee Mo with the DFT curve 
for bec Mo obtained using exactly the same procedures; we also show the fee melting curve 
of Belonoshko et al. 7 obtained using the Z method. 3 ^ We see that the free energy corrections 
cause a downward shift of the melting curve for both bec and fee, but the shift is considerably 
greater for fee. Consequently, the fee melting curve lies below the bec curve. In Table II, we 
show the values of terms (AU) l T s e{ and ((SAU) 2 ) re f (see Eq. (ITT]) ) which are required for the 
calculation of the free energy differences in question. The finding that the fee melting curve 
lies below the bec melting curve means that the free energy of fee must be higher than that 
of bec in the high-T region just below the melting curves. This confirms the conclusions 

23 



from our Bain-path and anharmonic calculations. However, our results are not consistent 
with those of Belonoshko et a/.—, whose Z-method calculations indicate that the fee melting 
curve lies above the bec melting curve. 

B. The Z method 

The electronic-structure methods used in our reference-coexistence calculations and in 
the Z-method calculations of Belonoshko et a/.- are essentially the same (PAW with the 
VASP code), so the contradictory conclusions about the relation between the bec and fee 
melting curves must originate in differences between the statistical-mechanical methods. 
The Z method has been validated by testing it against known results for the Lennard- Jones 
and other systems, using MD simulations on large systems of up to 32, 000 atoms with 
simulation times of ~ 60 ps^ 6 -. However, the DFT Z-method simulations of Belonoshko et 
al. 7 on the melting of Mo employed much smaller systems (from 32 to 108 atoms for fee, and 
from 54 to 128 atoms for bec), and very short simulation times of ~ 3 ps. It is therefore a 
natural question whether the use of such short simulations on such small systems might be 
the cause of the discrepancy. We have very recently investigated in detail the dependence 
of Z-method errors on system size and simulation time, and our findings shed light on this 
question^. Guided by this, we have performed our own DFT Z-method calculations on the 
melting of Mo, and we report the results here. 

The Z method is based on the phenomenon of homogeneous melting of a superheated 
solid 36 . The idea is that if an MD simulation is performed at constant total energy E and 
volume V (microcanonical ensemble) starting from the perfect crystal (all atoms on regular- 
lattice sites), then after the solid has thermally equilibrated at some temperature T so \ it will 
subsequently melt only if T sol exceeds a superheating limit T LS . Evidence was presented in 
Ref. I36J that, as the temperature T sol tends to T LS from above, the temperature Tn q and 
pressure Pu q of the liquid formed by homogeneous melting tend to a point on the melting 
curve. Our recent investigation of homogeneous melting^ focused on the waiting time r w , 
i.e. the time that elapses before the initial solid at temperature T so \ > Tls melts. In order to 
gather statistics about r w , for each system size (number of atoms N) with specified density 
N/V, and for each value of total energy (equivalently, for each equilibrated solid temperature 
^soi), we performed several hundred statistically independent simulations differing only in 
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T™ { (K) 


(AU)l s pf /N (eV/atom) 

\ / rci f \ / ' 


y((5AU) 2 } rci /N (eV/atom) 

J, < \ \ / / 1 , -'- L / \ / / 

Solid Liquid 


T£ l (K) 

III V / 


3200 


-0.057(2) 


0.024(2) 0.032(2) 


2249 


6325 


-0.108(2) 


0.044(2) 0.041(2) 


4910 


7625 


-0.070(2) 


0.015(2) 0.027(2) 


6690 



TABLE II: Difference (AU)[ s ei = (AU) l rci - (AU) s Ic{ between the liquid and fee solid thermal 
averages of the difference AU = Uai — U re f of ab initio and reference energies, and thermal averages 
in solid and liquid ((5AU) 2 ) TC { of the squared fluctuations of 5AU = AU — (Af7) re f, with averages 
evaluated in the reference system and normalized by dividing by the number of atoms N. Melting 
temperatures for the reference and ab initio systems are also reported. 
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the random velocities assigned at the start of the simulation. The key conclusions were 
that (a) r w is a stochastic quantity having a roughly exponential probability distribution; 
(b) its mean value (r w ) lengthens rapidly as T so \ — > T LS , being roughly proportional to 
l/(^soi — 2~ls) 2 ; ( c ) ( r w) also increases as the size of the system decreases, the dependence 
being roughly 1/N. We noted that if the total simulation time t S i m is much shorter than 
(t w ), then melting is unlikely to be observed even when T so \ > Tls- This means that if 
T m is estimated by performing simulations of fixed length t sim and seeking the lowest T so \ 
and Tii q for which melting is observed, then Tu q will inevitably overestimate T m , and the 
overestimation will become worse as the system size is reduced. As an indication of the 
difficulties, the results of Ref. |59| suggest that, for a transition metal with a system size of 
N = 100 and a simulation time t sim = 3 ps (values similar to those used by Belonoshko et 
al. 7 ), the overestimation could well be ~ 2000 K. Since the overestimation may differ for 
different crystal structures, it is clear that the Z method cannot be used to compare the 
melting temperatures of different crystal phases unless large enough systems are simulated 
for long enough times. 

To illustrate this point, we have performed our own DFT Z-method simulations on bcc 
and fee Mo, using systems of 250 atoms for bcc and 256 atoms for fee and simulation times 
of at least 12 ps (these are greater than the typical values used in Ref. Q by factors of 
2.5 and 4 respectively). The values of the final T and P in our simulations are reported 
in Fig. [TY] The results indicate that the fee crystal melts at a lower T than bcc, so that 
fee is thermodynamically less stable than bcc, as expected from our reference-coexistence 
calculations and from the free energies from our Bain-path and anharmonic calculations. 
We note that the conclusions from the present Z-method simulations are the opposite of the 



Z-method results of Ref. 



7|. This supports the suggestion that the earlier Z-method work 



employed simulations that were too short on systems that were too small. 



VI. DISCUSSION AND CONCLUSIONS 



Our results suggest that the fee and hep structures cannot be stable high-T phases of Mo 
in the pressure range < P < 600 GPa. We have shown that they would be more stable 
than bcc in the range 350 < P < 600 GPa and T > 5000 K in the harmonic approximation, 
as already found by Belonoshko et al.— and Zeng et al.— However, we find that anharmonic 
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330 335 340 345 350 

P (GPa) 

FIG. 11: Estimation of melting temperatures of bcc and fee Mo in the pressure region P ~ 345 
GPa. Round black points and red square points show final average (P,T) values from constant 
energy MD simulations (250 atoms for bcc, 256 atoms for fee) starting from the perfect lattice. On 
the left hand branches, melting does not occur within the duration of the simulations (at least 12 
ps, see text); on the right hand branches melting has occurred, and the (P,T) values refer to the 
liquid. 

contributions to the free energy substantially change the picture. Our most direct evidence 
for this in the case of fee comes from thermodynamic integration along the Bain path, which 
indicates that fee is thermodynamically less stable than bcc in the region where harmonic 
theory predicts the opposite; furthermore, fee appears to be elastically unstable at high T 
and P < 300 GPa. The Bain-path approach has the attractive feature that it relies only 
on completely standard first-principles MD, and the calculations are easily repeatable by 
other researchers. The existence of large anharmonic contributions, which crucially change 
the high-T stability of fee and hep relative to bcc, is confirmed by our explicit calculation 
of these contributions. Further confirmation that fee is thermodynamically less stable than 
bcc at high T comes from our comparison of the fee and bcc melting curves. 

At first sight, it might seem unexpected that anharmonicity stabilises bcc more than 
fee and hep. After all, fee and hep are the structures that go harmonically unstable at 
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P < 350 GPa, and intuition might suggest that below this pressure there could be large, 
anharmonically stabilised vibrations, which would have a large entropy. However, we have 
seen that the fee structure at high T is not vibrationally unstable, at least with the sizes 
of simulation cell that we have used, so presumably phonons that would be harmonically 
unstable are stiffened by anharmonic effects, so that their entropy is actually reduced. In 
fact, Asker et al. have shown recently that electronic thermal excitations have the effect 
of increasing the phonon frequencies of fee Mo (see Fig. 2 in Ref. 5JJ). Furthermore, 
electronic thermal excitations appear also to further stabilize the bec structure over fee. As 
we know from previous work, the general effect of high T is to smooth the peaks and valleys 
of the zero-temperature electronic DOS. However, in the bec structure the population of 
electronic states on the region near the Fermi energy is enhanced while in the fee structure 
it is depleted. This has the overall effect of enhancing the electronic entropy of the bec 
structure with respect to that of fee. A similar argument has already been suggested by 
Asker et a/.— for explaining the stability properties of Mo at low P and high T. It is also 
worth mentioning that in a recent study where we have developed a tight-binding model 
for Mo based on DFT data and used it to calculate anharmonic free energies over wide 
P — T intervals, no stabilization of the fee structure over bec is observed.— The effect of 
anharmonicity on the thermodynamic functions of the closely analogous element W has been 
discussed recently by Ozolins.— 

Our finding that the fee melting curve is below the bec curve, supported by our own 
Z-method calculations, is not consistent with the Mo melting curves deduced by Belonoshko 
et al— from their Z-method work. We have noted that one of the difficulties faced by the 
Z method concerns time scales. When the temperature T so \ of the initially thermalised 
crystal exceeds the superheating limit Tls, then in constant-energy MD the system will 
eventually melt, but the waiting time r w before this occurs may be tens of ps or even more 
if T so i is near Tls, so that l° n g simulations are needed if the method is to be reliable.— 
The time-scale problem appears to become worse for small systems. The evidence we have 
presented indicates that the simulation times of only ~ 3 ps used in the earlier Z-method 
work- were too short to yield reliable results. The longer simulations of at least 12 ps that 
we use here should give better results, but even so the bec melting temperature that we 
obtain is a significant overestimate compared with the values from our reference-coexistence 
calculations. It would clearly be desirable to repeat the Z-method calculations with still 
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longer runs on larger systems. However, the present simulations do serve the useful purpose 
of showing that the Z-method predictions for the relative melting temperatures of bcc and 
fee Mo can be consistent with our much more extensive and detailed results from free-energy 
calculations. 

The very recent DAC work of Dewaele et al.— on the melting of Ta makes it clear that 
very careful attention must be paid to experimental procedures if reliable results are to 
be obtained for high-P/high-T phase boundaries, and we believe that a cautious attitude 
should be adopted towards the existing DAC evidence^ for low melting curves in high- 
P Mo. Nevertheless, the shock data on Mo seem to require a crystallographic boundary 
somewhere in the region where quasiharmonic calculations indicate a transition from bcc 
to fee or hep. When we began the present work, we did not expect that the inclusion of 
anharmonicity would cause the bcc-fcc and bec-hep boundaries to disappear. Because we 
were initially sceptical of our findings, we felt it essential to confirm them in the ways that 
we have described. Our current belief is that efforts should be continued to search for other 
candidate crystal structures which might be thermodynamically more stable than bcc in the 
high-P/high-T region. 

In conclusion, our results suggest that the high-P/highT solid phase of Mo indicated 
by shock experiments is not fee or hep, but we do not rule out the possibility of other 
stable high-T crystal phases. Our results also suggest that the use of the quasiharmonic 
approximation should not be uncritically accepted in the first-principles search for other 
candidate crystal structures. 

Acknowledgments 

The work was supported by EPSRC-GB Grant No. EP/C534360, which was 50% funded 
by DSTL(MOD). The work was conducted as part of a EURYI scheme award to DA as 
provided by EPSRC-GB (see www.esf.org/euryi). The authors acknowledge an allocation 
time on the HECToR Supercomputer UK Facilities, some of which was provided by the 
UKCP consortium. We also benefited from the UCL Legion High Performance Computing 
Facilities and associated support services in the completion of this work. The authors thank 



29 



Prof. G. Kresse for help with PAW and Prof. C. J. Pickard for helpful discussions. 



1 C. Cazorla, M. J. Gillan, S. Taioli and D. Alfe, J. Chem. Phys. 126, 194502 (2007) 

2 S. Taioli, C. Cazorla, M. J. Gillan and D. Alfe, Phys. Rev. B 75, 214103 (2007) 

3 D. Errandonea, B. Schwager, R. Ditz, C. Gessman, R. Boehler and M. Ross, Phys. Rev. B 63, 
132104 (2001) 

4 D. Alfe, G. D. Price and M. J. Gillan, Phys. Rev. B 65, 165118 (2002) 

5 D. Errandonea, Physica B 357, 356 (2005) 

6 D. Santamaria-Perez, M. Ross, D. Errandonea, G. D. Mukherjee, M. Mezouar and R. Boehler, 
J. Chem. Phys. 130, 124509 (2009) 

7 A. B. Belonoshko, L. Burakovsky, S. P. Chen, B. Johansson, A. S. Mikhaylushkin, D. L. Preston, 
S. I. Simak and D. C. Swift, Phys. Rev. Lett. 100, 135701 (2008); Phys. Rev. Lett. 101, 049602 
(2008) 

8 T. S. Duffy, Rep. Prog. Phys., 68, 1811 (2005) 

9 L. Burakovsky, S. P. Chen, D. L. Preston, A. B. Belonoshko, A. Rosengren, A. S. Mikhaylushkin, 
S. I. Simak and J. A. Moriarty, Phys. Rev. Lett., 104, 255702 (2010) 

A. Dewaele, M. Mezouar, N. Guignot and P. Loubeyre, Phys. Rev. Lett., 104, 255701 (2010) 

1 C. J. Wu, P. Soderlind, J. N. Glosli and J. E. Klepeis, Nature Materials, 8, 223 (2009) 

2 Z.-L Liu, L.-C. Cai, X.-R. Chen and F.-Q. Jing, Phys. Rev. B, 77, 024103 (2008) 

3 R. S. Hixson and J. N. Fritz, J. Appl. Phys. 71, 1721 (1992) 

4 A. C. Mitchell and W. J. Nellis, J. Appl. Phys. 52, 3363 (1981) 

5 C. Cazorla, D. Alfe and M. J. Gillan , Phys. Rev. Lett. 101, 049601 (2008) 

6 Z.-Y. Zeng, C.-E. Hu, L.-C. Cai, X.-R. Chen and F.-Q. Jing, J. Phys. Chem. B, 114, 298 (2010) 

7 S. Japel, B. Schwager, R. Boehler and M. Ross, Phys. Rev. Lett. 95, 167801 (2005) 

8 D. Errandonea, M. Somayazulu, D. Hausermann and H. K. Mao, J. Phys.: Condens. Matter 
15, 7635 (2003) 

9 A. B. Belonoshko, S. I. Simak, A. E. Kochetov, B. Johansson, L. Burakovsky and D. L. Preston, 
Phys. Rev. Lett. 92, 195701 (2004) 

;0 L. Stixrude, R. E. Cohen and D. J. Singh, Phys. Rev. B 50, 6442 (1994) 

;1 L. Vocadlo, J. Brodholt, D. Alfe, M. J. Gillan and G. D. Price, Phys. of Earth and Planet. Int. 

30 



117, 123 (2000) 

22 S. Taioli, C. Cazorla, M. J. Gillan and D. Alfe, J. Phys.:Conf. Series 121, 012010 (2008) 

23 C. Cazorla, M. J. Gillan, S. Taioli and D. Alfe, J. Phys.:Conf. Series 121, 012009 (2008) 

24 D. Alfe, G. Kresse and M. J. Gillan, Phys. Rev. B, 61, 132 (2000) 

25 L. Vocadlo, I. G. Wood, M. J. Gillan, J. Brodholt, D. P. Dobson, G. D. Price and D. Alfe, Phys. 
Earth Planet. Int. 170, 52 (2008) 

26 H. K. Mao et al, Science 292, 914 (2001) 

27 C. Cazorla, D. Alfe and M. J. Gillan, Phys. Rev. B 77, 224103 (2008) 

28 G. A. de Wijs, G. Kresse and M. J. Gillan, Phys. Rev. B 57, 8223 (1998) 

29 D. Alfe, M. J. Gillan and G. D. Price, Nature 401, 462 (1999) 

30 A. Dewaele, M. Mezouar, N. Guignot and P. Loubeyre, Phys. Rev. B 76, 144106 (2007) 

31 L. Vocadlo and D. Alfe, Phys. Rev. B 65, 214105 (2002) 

32 O. Sugino and R. Car, Phys. Rev. Lett. 74, 1823 (1995) 

33 M. J. Gillan, D. Alfe, J. Brodholt, L. Vocadlo and G. D. Price, Rep. Prog. Phys. 69, 2365 (2006) 

34 L. Vocadlo, D. Alfe, M. J. Gillan and G. D. Price, J. Chem. Phys. 120, 2872 (2004) 

35 S. M. Davis, A. B. Belonoshko, B. Johansson, N. V. Skorodumova, and A. C. T. van Duin, J. 
Chem. Phys. 129, 194508 (2008) 

36 A. B. Belonoshko, N. V. Skorodumova, A. Rosengren and B. Johansson, Phys. Rev. B 73, 
012201 (2006) 

37 O. H. Nielsen and R. M. Martin, Phys. Rev. B 32, 3792 (1985) 

38 D. Alfe and M. J. Gillan, Phys. Rev. Lett. 81, 5161 (1998) 

39 L. Vocadlo, D. Alfe, M. J. Gillan, I. G. Wood, J. P. Brodholt and G. D. Price, Nature 424, 536 
(2003) 

40 D. Alfe, L. Vocadlo, M. J. Gillan and G. D. Price, J. Phys.: Condens. Matt. 16, S973 (2004) 

41 D. Alfe, M. J. Gillan and G. D. Price, J. Chem. Phys. 116, 6170 (2002) 

42 P. E. Blochl, Phys. Rev. B 50, 17953 (1994) 

43 G. Kresse and J. Furthmuller, Phys. Rev. B 54, 11169 (1996) 

44 J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996) 

45 H. J. Monkhorst and J. D. Pack, Phys. Rev. B 13, 5188 (1976) 

46 N. D. Mermin, Phys. Rev.137, A1441 (1965) 

47 D. Alfe, G. D. Price and M. J. Gillan, Phys. Rev. B 64, 045123 (2001) 

31 



48 G. Kresse, J. Furthmiiller and J. Hafner, Europhys. Lett. 32, 729 (1995) 

49 D. Alfe, Program available at http:/ /chianti.geol.ucl.ac.uk/^dario (1998); D. Alfe, Comp. Phys. 



Comm. 180, 2622 (2009) 

50 F. Birch, J. Geophys. Res. 83, 1257 (1978) 

51 C. Asker, A. B. Belonoshko, A. S. Mikhaylushkin and I. A. Abrikosov, Phys. Rev. B 77, 
220102(R) (2008) 

52 We note here the caveat that should be borne in mind for all types of thermodynamic integration: 
caution should be exercised if the path of integration passes through regions of thermodynam- 
ically unstable states, or if the thermally averaged integrand (stress, in the present case) is a 
discontinuous function of the integration variable (strain). The results we present show that 
the thermally averaged stress tensor varies continuously along the Bain path, so that there is 
no indication of instabilities on the path. Furthermore, the correctness of Bain-path integration 
will be confirmed a posteriori by the fact that the free-energy differences deduced from it are 
corroborated by explicit calculation of the anharmonic free energies (see Sec. lIVp . 

53 M. S. Daw and M. I. Baskes, Phys. Rev. B 29, 6443 (1984) 

54 M. W. Finnis and J. E. Sinclair, Phil. Mag. A 50, 45 (1984) 

55 We have also checked the effect on the fcc-bcc free energy difference when we reduce the PAW 
core radius. We reduced the core radius from 2.7 a.u. to 1.6 a.u. and we used first-order perturba- 
tion theory to estimate the changes of free energies of fee and bec at the state V = 8.26 A 3 /atom 
and T = 9000 K. We found that these free energy changes amount to only ~ 10 meV/atom, 
which is much smaller than our computed free energy differences between fee and bcc. 

56 We note that the use of harmonic reference systems in calculating the free energy of anharmonic 
systems is very well established. A few of the many examples are: D. Frenkel and A. J. C. Ladd, 
J. Chem. Phys. 81, 3188 (1984); S. M. Foiles, Phys. Rev. B 49, 14930 (1994); O. Sugino and 
R. Car, Phys. Rev. Lett. 74, 1823 (1995); U. Hansen, P. Vogl and V. Fiorentini, Phys. Rev. B 
60, 5055 (1999). 

57 D. Alfe, Phys. Rev. B 79, 060101(R) (2009) 

58 E. Sola and D. Alfe, Phys. Rev. Lett. 103, 078501 (2009) 

59 D. Alfe, C. Cazorla and M. J. Gillan, J. Chem. Phys., submitted, arXiv:1104.2147l 

60 V. Ozolins, Phys. Rev. Lett. 102, 065702 (2009) 



61 This point was emphasised by the originators of the Z-method, noting in Ref. 

32 



36 that "...in this 



method both a rather large number of atoms and long runs are needed...". 
C. Cazorla, D. Alfe and M. J. Gillan, J. Chem. Phys. 130, 174707 (2009) 



33 



